<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>slasq2.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SLASQ2.1"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>( N, Z, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Modified to call <a name="SLAZQ3.7"></a><a href="slazq3.f.html#SLAZQ3.1">SLAZQ3</a> in place of <a name="SLASQ3.7"></a><a href="slasq3.f.html#SLASQ3.1">SLASQ3</a>, 13 Feb 03, SJH.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            INFO, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               Z( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SLASQ2.19"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a> computes all the eigenvalues of the symmetric positive 
</span><span class="comment">*</span><span class="comment">  definite tridiagonal matrix associated with the qd array Z to high
</span><span class="comment">*</span><span class="comment">  relative accuracy are computed to high relative accuracy, in the
</span><span class="comment">*</span><span class="comment">  absence of denormalization, underflow and overflow.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  To see the relation of Z to the tridiagonal matrix, let L be a
</span><span class="comment">*</span><span class="comment">  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
</span><span class="comment">*</span><span class="comment">  let U be an upper bidiagonal matrix with 1's above and diagonal
</span><span class="comment">*</span><span class="comment">  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
</span><span class="comment">*</span><span class="comment">  symmetric tridiagonal to which it is similar.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Note : <a name="SLASQ2.30"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a> defines a logical variable, IEEE, which is true
</span><span class="comment">*</span><span class="comment">  on machines which follow ieee-754 floating-point standard in their
</span><span class="comment">*</span><span class="comment">  handling of infinities and NaNs, and false otherwise. This variable
</span><span class="comment">*</span><span class="comment">  is passed to <a name="SLAZQ3.33"></a><a href="slazq3.f.html#SLAZQ3.1">SLAZQ3</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N     (input) INTEGER
</span><span class="comment">*</span><span class="comment">        The number of rows and columns in the matrix. N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Z     (workspace) REAL array, dimension (4*N)
</span><span class="comment">*</span><span class="comment">        On entry Z holds the qd array. On exit, entries 1 to N hold
</span><span class="comment">*</span><span class="comment">        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
</span><span class="comment">*</span><span class="comment">        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
</span><span class="comment">*</span><span class="comment">        N &gt; 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
</span><span class="comment">*</span><span class="comment">        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
</span><span class="comment">*</span><span class="comment">        shifts that failed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO  (output) INTEGER
</span><span class="comment">*</span><span class="comment">        = 0: successful exit
</span><span class="comment">*</span><span class="comment">        &lt; 0: if the i-th argument is a scalar and had an illegal
</span><span class="comment">*</span><span class="comment">             value, then INFO = -i, if the i-th argument is an
</span><span class="comment">*</span><span class="comment">             array and the j-entry had an illegal value, then
</span><span class="comment">*</span><span class="comment">             INFO = -(i*100+j)
</span><span class="comment">*</span><span class="comment">        &gt; 0: the algorithm failed
</span><span class="comment">*</span><span class="comment">              = 1, a split was marked by a positive value in E
</span><span class="comment">*</span><span class="comment">              = 2, current block of Z not diagonalized after 30*N
</span><span class="comment">*</span><span class="comment">                   iterations (in inner while loop)
</span><span class="comment">*</span><span class="comment">              = 3, termination criterion of outer while loop not met 
</span><span class="comment">*</span><span class="comment">                   (program created more than N unreduced blocks)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">  Local Variables: I0:N0 defines a current unreduced segment of Z.
</span><span class="comment">*</span><span class="comment">  The shifts are accumulated in SIGMA. Iteration count is in ITER.
</span><span class="comment">*</span><span class="comment">  Ping-pong is controlled by PP (alternates between 0 and 1).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               CBIAS
      PARAMETER          ( CBIAS = 1.50E0 )
      REAL               ZERO, HALF, ONE, TWO, FOUR, HUNDRD
      PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
     $                     TWO = 2.0E0, FOUR = 4.0E0, HUNDRD = 100.0E0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            IEEE
      INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, 
     $                   N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
      REAL               D, DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, E,
     $                   EMAX, EMIN, EPS, OLDEMN, QMAX, QMIN, S, SAFMIN,
     $                   SIGMA, T, TAU, TEMP, TOL, TOL2, TRACE, ZMAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLAZQ3.86"></a><a href="slazq3.f.html#SLAZQ3.1">SLAZQ3</a>, <a name="SLASRT.86"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>, <a name="XERBLA.86"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            <a name="ILAENV.89"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      REAL               <a name="SLAMCH.90"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="ILAENV.91"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>, <a name="SLAMCH.91"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN, REAL, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">      
</span><span class="comment">*</span><span class="comment">     Test the input arguments.
</span><span class="comment">*</span><span class="comment">     (in case <a name="SLASQ2.99"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a> is not called by <a name="SLASQ1.99"></a><a href="slasq1.f.html#SLASQ1.1">SLASQ1</a>)
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      EPS = <a name="SLAMCH.102"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Precision'</span> )
      SAFMIN = <a name="SLAMCH.103"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
      TOL = EPS*HUNDRD
      TOL2 = TOL**2
<span class="comment">*</span><span class="comment">
</span>      IF( N.LT.0 ) THEN
         INFO = -1
         CALL <a name="XERBLA.109"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLASQ2.109"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, 1 )
         RETURN
      ELSE IF( N.EQ.0 ) THEN
         RETURN
      ELSE IF( N.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        1-by-1 case.
</span><span class="comment">*</span><span class="comment">
</span>         IF( Z( 1 ).LT.ZERO ) THEN
            INFO = -201
            CALL <a name="XERBLA.119"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLASQ2.119"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, 2 )
         END IF
         RETURN
      ELSE IF( N.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        2-by-2 case.
</span><span class="comment">*</span><span class="comment">
</span>         IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
            INFO = -2
            CALL <a name="XERBLA.128"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLASQ2.128"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, 2 )
            RETURN
         ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
            D = Z( 3 )
            Z( 3 ) = Z( 1 )
            Z( 1 ) = D
         END IF
         Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
         IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
            T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) 
            S = Z( 3 )*( Z( 2 ) / T )
            IF( S.LE.T ) THEN
               S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
            ELSE
               S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
            END IF
            T = Z( 1 ) + ( S+Z( 2 ) )
            Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
            Z( 1 ) = T
         END IF
         Z( 2 ) = Z( 3 )
         Z( 6 ) = Z( 2 ) + Z( 1 )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for negative data and compute sums of q's and e's.
</span><span class="comment">*</span><span class="comment">
</span>      Z( 2*N ) = ZERO
      EMIN = Z( 2 )
      QMAX = ZERO
      ZMAX = ZERO
      D = ZERO
      E = ZERO
<span class="comment">*</span><span class="comment">
</span>      DO 10 K = 1, 2*( N-1 ), 2
         IF( Z( K ).LT.ZERO ) THEN
            INFO = -( 200+K )
            CALL <a name="XERBLA.165"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLASQ2.165"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, 2 )
            RETURN
         ELSE IF( Z( K+1 ).LT.ZERO ) THEN
            INFO = -( 200+K+1 )
            CALL <a name="XERBLA.169"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLASQ2.169"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, 2 )
            RETURN
         END IF
         D = D + Z( K )
         E = E + Z( K+1 )
         QMAX = MAX( QMAX, Z( K ) )
         EMIN = MIN( EMIN, Z( K+1 ) )
         ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
   10 CONTINUE
      IF( Z( 2*N-1 ).LT.ZERO ) THEN
         INFO = -( 200+2*N-1 )
         CALL <a name="XERBLA.180"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLASQ2.180"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, 2 )
         RETURN
      END IF
      D = D + Z( 2*N-1 )
      QMAX = MAX( QMAX, Z( 2*N-1 ) )
      ZMAX = MAX( QMAX, ZMAX )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for diagonality.
</span><span class="comment">*</span><span class="comment">
</span>      IF( E.EQ.ZERO ) THEN
         DO 20 K = 2, N
            Z( K ) = Z( 2*K-1 )
   20    CONTINUE
         CALL <a name="SLASRT.193"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'D'</span>, N, Z, IINFO )
         Z( 2*N-1 ) = D
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      TRACE = D + E
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for zero data.
</span><span class="comment">*</span><span class="comment">
</span>      IF( TRACE.EQ.ZERO ) THEN
         Z( 2*N-1 ) = ZERO
         RETURN
      END IF
<span class="comment">*</span><span class="comment">         
</span><span class="comment">*</span><span class="comment">     Check whether the machine is IEEE conformable.
</span><span class="comment">*</span><span class="comment">         
</span>      IEEE = <a name="ILAENV.209"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 10, <span class="string">'<a name="SLASQ2.209"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, <span class="string">'N'</span>, 1, 2, 3, 4 ).EQ.1 .AND.
     $       <a name="ILAENV.210"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 11, <span class="string">'<a name="SLASQ2.210"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>'</span>, <span class="string">'N'</span>, 1, 2, 3, 4 ).EQ.1      
<span class="comment">*</span><span class="comment">         
</span><span class="comment">*</span><span class="comment">     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
</span><span class="comment">*</span><span class="comment">
</span>      DO 30 K = 2*N, 2, -2
         Z( 2*K ) = ZERO 
         Z( 2*K-1 ) = Z( K ) 
         Z( 2*K-2 ) = ZERO 
         Z( 2*K-3 ) = Z( K-1 ) 
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      I0 = 1
      N0 = N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Reverse the qd-array, if warranted.
</span><span class="comment">*</span><span class="comment">
</span>      IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
         IPN4 = 4*( I0+N0 )
         DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
            TEMP = Z( I4-3 )
            Z( I4-3 ) = Z( IPN4-I4-3 )
            Z( IPN4-I4-3 ) = TEMP
            TEMP = Z( I4-1 )
            Z( I4-1 ) = Z( IPN4-I4-5 )
            Z( IPN4-I4-5 ) = TEMP
   40    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initial split checking via dqd and Li's test.
</span><span class="comment">*</span><span class="comment">
</span>      PP = 0
<span class="comment">*</span><span class="comment">
</span>      DO 80 K = 1, 2
<span class="comment">*</span><span class="comment">
</span>         D = Z( 4*N0+PP-3 )
         DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
            IF( Z( I4-1 ).LE.TOL2*D ) THEN
               Z( I4-1 ) = -ZERO
               D = Z( I4-3 )
            ELSE
               D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
            END IF
   50    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        dqd maps Z to ZZ plus Li's test.
</span><span class="comment">*</span><span class="comment">
</span>         EMIN = Z( 4*I0+PP+1 )
         D = Z( 4*I0+PP-3 )
         DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
            Z( I4-2*PP-2 ) = D + Z( I4-1 )
            IF( Z( I4-1 ).LE.TOL2*D ) THEN
               Z( I4-1 ) = -ZERO
               Z( I4-2*PP-2 ) = D
               Z( I4-2*PP ) = ZERO
               D = Z( I4+1 )
            ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
     $               SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
               TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
               Z( I4-2*PP ) = Z( I4-1 )*TEMP
               D = D*TEMP
            ELSE
               Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
               D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
            END IF
            EMIN = MIN( EMIN, Z( I4-2*PP ) )
   60    CONTINUE 
         Z( 4*N0-PP-2 ) = D
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now find qmax.
</span><span class="comment">*</span><span class="comment">
</span>         QMAX = Z( 4*I0-PP-2 )
         DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
            QMAX = MAX( QMAX, Z( I4 ) )
   70    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Prepare for the next iteration on K.
</span><span class="comment">*</span><span class="comment">
</span>         PP = 1 - PP
   80 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialise variables to pass to <a name="SLAZQ3.290"></a><a href="slazq3.f.html#SLAZQ3.1">SLAZQ3</a>
</span><span class="comment">*</span><span class="comment">
</span>      TTYPE = 0
      DMIN1 = ZERO
      DMIN2 = ZERO
      DN    = ZERO
      DN1   = ZERO
      DN2   = ZERO
      TAU   = ZERO
<span class="comment">*</span><span class="comment">
</span>      ITER = 2
      NFAIL = 0
      NDIV = 2*( N0-I0 )
<span class="comment">*</span><span class="comment">
</span>      DO 140 IWHILA = 1, N + 1
         IF( N0.LT.1 ) 
     $      GO TO 150
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        While array unfinished do 
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        E(N0) holds the value of SIGMA when submatrix in I0:N0
</span><span class="comment">*</span><span class="comment">        splits from the rest of the array, but is negated.
</span><span class="comment">*</span><span class="comment">      
</span>         DESIG = ZERO
         IF( N0.EQ.N ) THEN
            SIGMA = ZERO
         ELSE
            SIGMA = -Z( 4*N0-1 )
         END IF
         IF( SIGMA.LT.ZERO ) THEN
            INFO = 1
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Find last unreduced submatrix's top index I0, find QMAX and
</span><span class="comment">*</span><span class="comment">        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
</span><span class="comment">*</span><span class="comment">
</span>         EMAX = ZERO 
         IF( N0.GT.I0 ) THEN
            EMIN = ABS( Z( 4*N0-5 ) )
         ELSE
            EMIN = ZERO
         END IF
         QMIN = Z( 4*N0-3 )
         QMAX = QMIN
         DO 90 I4 = 4*N0, 8, -4
            IF( Z( I4-5 ).LE.ZERO )
     $         GO TO 100
            IF( QMIN.GE.FOUR*EMAX ) THEN
               QMIN = MIN( QMIN, Z( I4-3 ) )
               EMAX = MAX( EMAX, Z( I4-5 ) )
            END IF
            QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
            EMIN = MIN( EMIN, Z( I4-5 ) )
   90    CONTINUE
         I4 = 4 
<span class="comment">*</span><span class="comment">
</span>  100    CONTINUE
         I0 = I4 / 4
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Store EMIN for passing to <a name="SLAZQ3.350"></a><a href="slazq3.f.html#SLAZQ3.1">SLAZQ3</a>.
</span><span class="comment">*</span><span class="comment">
</span>         Z( 4*N0-1 ) = EMIN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Put -(initial shift) into DMIN.
</span><span class="comment">*</span><span class="comment">
</span>         DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong.
</span><span class="comment">*</span><span class="comment">
</span>         PP = 0 
<span class="comment">*</span><span class="comment">
</span>         NBIG = 30*( N0-I0+1 )
         DO 120 IWHILB = 1, NBIG
            IF( I0.GT.N0 ) 
     $         GO TO 130
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           While submatrix unfinished take a good dqds step.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLAZQ3.369"></a><a href="slazq3.f.html#SLAZQ3.1">SLAZQ3</a>( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
     $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
     $                   DN2, TAU )
<span class="comment">*</span><span class="comment">
</span>            PP = 1 - PP
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           When EMIN is very small check for splits.
</span><span class="comment">*</span><span class="comment">
</span>            IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
               IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
     $             Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
                  SPLT = I0 - 1
                  QMAX = Z( 4*I0-3 )
                  EMIN = Z( 4*I0-1 )
                  OLDEMN = Z( 4*I0 )
                  DO 110 I4 = 4*I0, 4*( N0-3 ), 4
                     IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
     $                   Z( I4-1 ).LE.TOL2*SIGMA ) THEN
                        Z( I4-1 ) = -SIGMA
                        SPLT = I4 / 4
                        QMAX = ZERO
                        EMIN = Z( I4+3 )
                        OLDEMN = Z( I4+4 )
                     ELSE
                        QMAX = MAX( QMAX, Z( I4+1 ) )
                        EMIN = MIN( EMIN, Z( I4-1 ) )
                        OLDEMN = MIN( OLDEMN, Z( I4 ) )
                     END IF
  110             CONTINUE
                  Z( 4*N0-1 ) = EMIN
                  Z( 4*N0 ) = OLDEMN
                  I0 = SPLT + 1
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         INFO = 2
         RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        end IWHILB
</span><span class="comment">*</span><span class="comment">
</span>  130    CONTINUE
<span class="comment">*</span><span class="comment">
</span>  140 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      INFO = 3
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     end IWHILA   
</span><span class="comment">*</span><span class="comment">
</span>  150 CONTINUE
<span class="comment">*</span><span class="comment">      
</span><span class="comment">*</span><span class="comment">     Move q's to the front.
</span><span class="comment">*</span><span class="comment">      
</span>      DO 160 K = 2, N
         Z( K ) = Z( 4*K-3 )
  160 CONTINUE
<span class="comment">*</span><span class="comment">      
</span><span class="comment">*</span><span class="comment">     Sort and compute sum of eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLASRT.430"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'D'</span>, N, Z, IINFO )
<span class="comment">*</span><span class="comment">
</span>      E = ZERO
      DO 170 K = N, 1, -1
         E = E + Z( K )
  170 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Store trace, sum(eigenvalues) and information on performance.
</span><span class="comment">*</span><span class="comment">
</span>      Z( 2*N+1 ) = TRACE 
      Z( 2*N+2 ) = E
      Z( 2*N+3 ) = REAL( ITER )
      Z( 2*N+4 ) = REAL( NDIV ) / REAL( N**2 )
      Z( 2*N+5 ) = HUNDRD*NFAIL / REAL( ITER )
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLASQ2.446"></a><a href="slasq2.f.html#SLASQ2.1">SLASQ2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
